home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / LIBIP / PMOM.C < prev    next >
C/C++ Source or Header  |  1999-09-11  |  10KB  |  346 lines
  1. /* 
  2.  * pmom.c P(oly)MOM(ents)
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. #include <stdio.h>
  10. #include <stdlib.h>
  11. #include <math.h>
  12. #include "ip.h"
  13.  
  14. #define    DENT_CUTOFF    0.001
  15. #define    TANGENT_LIMIT    1.0e06    /* indicates vertical axis        */
  16. #define    MU11_EPS    5.0e03         /* indicates alignment with axes */
  17.  
  18. #define ZERO        0.0
  19. #define FABS(a)        ((a) > ZERO ? (a) : - (a))
  20. #define    SQR(a)        ( (a)*(a) )
  21.  
  22. #define ON        1
  23. #define OFF        0
  24. #define SAVE_MOMENTS    ON
  25.  
  26. #define    DEBUG
  27. #define    CEN_DEBUG
  28. #define    SHOW_VERT
  29.  
  30. /*
  31.  * poly_moments()
  32.  *   DESCRIPTION:
  33.  *       poly_moments computes centroid, central moments, principal axes and
  34.  *       first two moment invariants
  35.  *     the centroid position is marked with a cross
  36.  *     the principle axis directions are plotted
  37.  *   ARGUMENTS:
  38.  *     nv(long) number of vertices in polygon
  39.  *     v(struct spoint *) pointer to struct containing the vertices
  40.  *     moments(float *) pointer to array of floats to hold calculated moments
  41.  *     imgIO(Image *) pointer to Image struct (see tiffimage.h)
  42.  *     value(int) pixel value to use for plotting
  43.  *   RETURN VALUE:
  44.  *     (struct spoint *) that contains the centroid value
  45.  */
  46.  
  47. struct spoint *
  48. poly_moments (long nv, struct spoint *v, float *moments, Image * imgIO, int value)
  49. {
  50.   int i, i_max, i_min;
  51.   int nvpp;
  52.   int size = 30;
  53.   int circle = 0, horizontal = 0, vertical = 0;
  54.  
  55.   double ximm, xi, yimm, yi;
  56.   double d_xy, d_min, d_max;
  57.  
  58.   double m00, m00_sum;
  59.   double m10, m10_sum, m01, m01_sum, m11, m11_sum;
  60.   double m20, m20_sum, m02, m02_sum;
  61.   double mu, mu11, mu20, mu02, dent;
  62.   double musq, phi_1, phi_2;
  63.   double rad_gyr;
  64.   double tg_tth, tg_th1, tg_th2, sq_root;
  65.   double mu02_div_mu20, mu11mu20_sign;
  66.   double x_major, y_major, x_minor, y_minor;
  67.   int quad_IV = OFF, quad_I = OFF;
  68.  
  69.   static struct spoint *pvc, vc;
  70.  
  71.   int color_index = 191;
  72.   pvc = &vc;
  73.  
  74.  
  75. #ifdef SHOW_VERT
  76.   printf ("\n...PMOM: vertex coordinates\n");
  77.   for (i = 0; i < (int) nv; i++)
  78.     printf ("...%d: (%3d, %3d)\n", i, (v + i)->x, (v + i)->y);
  79. #endif
  80.  
  81. /*
  82.  * compute moments
  83.  * refs:
  84.  *   Hu, IRE Trans. Inf. Theory, IT-8, 179-187 (1962)
  85.  *   Seul, Sammon and Monar, Rev. Sci. Instr. 62, 784-792 (1991)
  86.  */
  87.   m00 = m10 = m01 = 0.0;
  88.   m11 = m20 = m02 = 0.0;
  89.   nvpp = (int) nv + 1;
  90.   for (i = 1; i < nvpp; i++) {
  91.     ximm = (double) (v + i - 1)->x;
  92.     xi = (double) (v + i)->x;
  93.     yimm = (double) (v + i - 1)->y;
  94.     yi = (double) (v + i)->y;
  95.  
  96.     m00_sum = 0.5 * (yi * ximm - xi * yimm);
  97.     m00 += m00_sum;
  98.  
  99.     m10_sum = 0.5 * (xi + ximm) * m00_sum;
  100.     m10_sum -= 0.5 * ((yi - yimm) * (xi * xi + xi * ximm + ximm * ximm) / 6.0);
  101.     m10 += m10_sum;
  102.  
  103.     m01_sum = 0.5 * (yi + yimm) * m00_sum;
  104.     m01_sum += 0.5 * ((xi - ximm) * (yi * yi + yi * yimm + yimm * yimm) / 6.0);
  105.     m01 += m01_sum;
  106.  
  107.     m11_sum = 0.5 * m00_sum;
  108.     m11 += m11_sum * (2.0 * xi * yi + ximm * yi + xi * yimm + 2.0 * ximm * yimm) / 6.0;
  109.  
  110.     m20_sum = m00_sum * (xi * xi + xi * ximm + ximm * ximm) / 3.0;
  111.     m20_sum -= 0.5 * (yi - yimm) * (xi * xi * xi + xi * xi * ximm + xi * ximm * ximm + ximm * ximm * ximm) / 6.0;
  112.     m20 += m20_sum;
  113.  
  114.     m02_sum = m00_sum * (yi * yi + yi * yimm + yimm * yimm) / 3.0;
  115.     m02_sum += 0.5 * (xi - ximm) * (yi * yi * yi + yi * yi * yimm + yi * yimm * yimm + yimm * yimm * yimm) / 6.0;
  116.     m02 += m02_sum;
  117.  
  118.   }
  119.   printf ("\nmoments:\n");
  120.   printf ("...m00 = %lf\n", m00);
  121.   printf ("...m10 = %lf  m01 = %lf\n", m10, m01);
  122.   printf ("...m11 = %lf  m20 = %lf  m02 = %lf\n", m11, m20, m02);
  123.  
  124.  
  125. /*
  126.  * correct sign of raw moments if necessary
  127.  * ( -->negative moments are obtained from Zahn's representation
  128.  *      of an 'inner' boundary, traversed in CCW sense)
  129.  */
  130.   if (SIGN (m00) < ZERO) {
  131.     m00 *= -1.0;
  132.     m10 *= -1.0;
  133.     m01 *= -1.0;
  134.     m11 *= -1.0;
  135.     m20 *= -1.0;
  136.     m02 *= -1.0;
  137.   }
  138.  
  139.  
  140. /*
  141.  * centroid 
  142.  */
  143.   pvc->x = (int) (m10 / m00);
  144.   pvc->y = (int) (m01 / m00);
  145. #ifdef CEN_DEBUG
  146.   printf ("\nPOLY_MOMENT: centroid pos: (%3d, %3d)\n", pvc->x, pvc->y);
  147. #endif
  148.  
  149. /*
  150.  * mark centroid position with a cross
  151.  */
  152.   draw_cross (pvc->x, pvc->y, size, imgIO, value);
  153.  
  154.  
  155. /*
  156.  * find curvature points closest to and farthest from the centroid
  157.  */
  158.   i_max = i_min = 0;
  159.   d_max = 0.0;
  160.   d_min = 1000.0;
  161.   for (i = 1; i < nvpp; i++) {
  162.     xi = (double) (v + i)->x;
  163.     yi = (double) (v + i)->y;
  164.     d_xy = sqrt (SQR (xi - pvc->x) + SQR (yi - pvc->y));
  165.     if (d_xy > d_max) {
  166.       d_max = d_xy;
  167.       i_max = i;
  168.     }
  169.     if (d_xy < d_min) {
  170.       d_min = d_xy;
  171.       i_min = i;
  172.     }
  173.   }
  174. #ifdef DEBUG
  175.   printf ("\n...d_min = %lf  d_max = %lf\n", d_min, d_max);
  176. #endif
  177.  
  178. /*
  179.  * central moments and radius of gyration
  180.  * ref: Gonzalez & Wintz, "Digital Image Processing", chapt. 8.3
  181.  *
  182.  * at the expense of efficiency, work with normalized raw moments
  183.  * where possible, to keep numbers small
  184.  */
  185.   mu = m00;
  186.   mu11 = m11 - m00 * (m10 / m00) * (m01 / m00);
  187.   mu20 = m20 - m00 * SQR (m10 / m00);
  188.   mu02 = m02 - m00 * SQR (m01 / m00);
  189.  
  190.   rad_gyr = sqrt (FABS (mu20 + mu02) / mu);
  191.   dent = FABS (((mu02 / mu20) - 1.0));
  192.  
  193.   printf ("\ncentral moments:\n");
  194.   printf ("...mu00 = %lf\n", mu);
  195.   printf ("...mu11 = %lf  mu20 = %lf  mu02 = %lf\n", mu11, mu20, mu02);
  196.   printf ("...radius of gyration: %lf\n", rad_gyr);
  197.   printf ("...(mu02 - mu20)/mu20 = %lf\n", dent);
  198.  
  199.   if ((dent < DENT_CUTOFF) && (FABS (mu11) < MU11_EPS)) {
  200.     circle = ON;
  201.     printf ("\n--->circle!!\n");
  202.   }
  203.  
  204. /*
  205.  * invariant moments
  206.  * (reference: Gonzalez & Wintz, "Digital Image Processing", chapt. 8.3)
  207.  */
  208.   musq = SQR (mu);
  209.   phi_1 = (mu20 + mu02) / musq;
  210.   phi_2 = SQR ((mu20 - mu02) / musq) + 4.0 * SQR (mu11 / musq);
  211.  
  212.   printf ("\ninvariant moments:\n");
  213.   printf ("...phi_1 = %lf  phi_2 = %lf\n", phi_1, phi_2);
  214.  
  215. /*
  216.  * determine principal axes
  217.  * ref: A. Rosenfeld, A. C. Kak, "Digital Picture Processing", vol.II 
  218.  */
  219.   mu02_div_mu20 = mu02 / mu20;
  220.   tg_tth = 2.0 * (mu11 / mu20) / (1.0 - mu02_div_mu20);
  221.   printf ("\n...tg_tth = %lf\n", tg_tth);
  222.  
  223. /*
  224.  * check limits
  225.  */
  226.   if (FABS (tg_tth) < 0.01) {
  227.     if (FABS (mu20) / FABS (mu02) >= 1.0)
  228.       horizontal = ON;
  229.     if (FABS (mu20) / FABS (mu02) <= 1.0)
  230.       vertical = ON;
  231.   }
  232.   if ((circle == ON) || (horizontal == ON) || (vertical == ON)) {
  233.     printf ("\n...circle = %d", circle);
  234.     printf ("...horizontal = %d", horizontal);
  235.     printf ("...vertical = %d", vertical);
  236.     printf ("\n-->principal axes along screen coordinate system\n");
  237.     if ((circle == ON) || (horizontal == ON)) {
  238.       x_major = d_min + (double) pvc->x;
  239.       y_major = (double) pvc->y;
  240.       x_minor = (double) pvc->x;
  241.       y_minor = d_min + (double) pvc->y;
  242.       tg_th1 = ZERO;
  243.       tg_th2 = TANGENT_LIMIT;
  244.     }
  245.     else if (vertical == ON) {
  246.       x_major = (double) pvc->x;
  247.       y_major = d_min + (double) pvc->y;
  248.       x_minor = d_min + (double) pvc->x;
  249.       y_minor = (double) pvc->y;
  250.       tg_th1 = TANGENT_LIMIT;
  251.       tg_th2 = ZERO;
  252.     }
  253.   }
  254.   else {
  255.     sq_root = sqrt (1.0 + 1.0 / SQR (tg_tth));
  256.  
  257. /*
  258.  * determine proper direction of principal (major) axis 
  259.  * (sign of tg_th1, below):
  260.  *   -->check whether mu11 has same or opposite sign as do mu20 and m02
  261.  */
  262.     if (SIGN (mu20) * SIGN (mu02) < ZERO) {
  263.       printf ("\n...mu20, mu02 differ in sign");
  264.       printf (" -->something wrong?\n");
  265.     }
  266.  
  267.     if ((mu11mu20_sign = SIGN (mu11) * SIGN (mu20)) > ZERO) {
  268.       quad_IV = ON;             /* quadr II->IV */
  269.  
  270.       if (mu02_div_mu20 <= 1.0)
  271.         tg_tth = -FABS (tg_tth);
  272.       else
  273.         tg_tth = FABS (tg_tth);
  274.  
  275.       tg_th1 = (-1.0 / tg_tth) - sq_root;
  276.       tg_th2 = (-1.0 / tg_tth) + sq_root;
  277.  
  278.     }
  279.     else if (mu11mu20_sign < ZERO) {
  280.       quad_I = ON;              /* quadr III->I */
  281.  
  282.       if (mu02_div_mu20 <= 1.0)
  283.         tg_tth = FABS (tg_tth);
  284.       else
  285.         tg_tth = -FABS (tg_tth);
  286.  
  287.       tg_th1 = (-1.0 / tg_tth) + sq_root;
  288.       tg_th2 = (-1.0 / tg_tth) - sq_root;
  289.  
  290.     }
  291.     else
  292.       printf ("\n...unknown condition encountered\n");
  293.  
  294.  
  295.     printf ("\nprincipal axes:\n");
  296.     printf ("...tg_tth = %lf  tg_th1 = %lf  tg_th2 = %lf\n",
  297.             tg_tth, tg_th1, tg_th2);
  298.  
  299. /*
  300.  * determine coordinates of endpoints of line segments (of length 
  301.  * d_min and 0.5*d_min) to be drawn in direction of principal axes
  302.  */
  303.     x_major = d_min / sqrt (1.0 + SQR (tg_th1));
  304.     x_minor = 0.5 * d_min / sqrt (1.0 + SQR (tg_th2));
  305.     if (quad_I == ON) {
  306.       y_major = (double) pvc->y - FABS (tg_th1) * x_major;
  307.       y_minor = (double) pvc->y + FABS (tg_th2) * x_minor;
  308.     }
  309.     if (quad_IV == ON) {
  310.       y_major = (double) pvc->y + FABS (tg_th1) * x_major;
  311.       y_minor = (double) pvc->y - FABS (tg_th2) * x_minor;
  312.     }
  313.     x_major += (double) pvc->x;
  314.     x_minor += (double) pvc->x;
  315.   }
  316.   printf ("\n...x_major = %lf  y_major = %lf\n", x_major, y_major);
  317.   printf ("...x_minor = %lf  y_minor = %lf\n", x_minor, y_minor);
  318.  
  319. /*
  320.  * mark principal axis directions
  321.  */
  322.   draw_line (pvc->x, pvc->y, (int) x_major, (int) y_major, imgIO, value);
  323.   draw_line (pvc->x, pvc->y, (int) x_minor, (int) y_minor, imgIO, value);
  324.  
  325. /*
  326.  * save paramters related to moments for later reference
  327.  */
  328.   *(moments + 0) = (float) SAVE_MOMENTS;
  329.   *(moments + 1) = (float) m00;
  330.   *(moments + 2) = (float) m10;
  331.   *(moments + 3) = (float) m01;
  332.   *(moments + 4) = (float) m11;
  333.   *(moments + 5) = (float) m20;
  334.   *(moments + 6) = (float) m02;
  335.   *(moments + 7) = (float) mu11;
  336.   *(moments + 8) = (float) mu20;
  337.   *(moments + 9) = (float) mu02;
  338.   *(moments + 10) = (float) dent;
  339.   *(moments + 11) = (float) phi_1;
  340.   *(moments + 12) = (float) phi_2;
  341.   *(moments + 13) = (float) tg_th1;
  342.   *(moments + 14) = (float) tg_th2;
  343.  
  344.   return (pvc);
  345. }
  346.